% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function h = plot_u(grid,base,U,varargin)
% h = plot_u(grid,base,U)
% h = plot_u(grid,base,U,k)
%
% Plot a scalar U (represented as a discontinuous field, i.e. as a two
% dimensional array). On each element, the solution is resampled on a
% local mesh, which is then plotted with my_pdesurf. If k is absent,
% the degree of the basis is used.

 % 0) preliminary checks
 if(nargin==3)
   k = base.k;
 else
   k = varargin{1};
 end
 
 % 1) resample the solution
 [Xi,Ti,Ui] = el_resample(grid,base,k,U);

 % 2) open the figure and plot the mesh
 h = figure;
 my_pdemesh(h,grid,'k');
 hold on

 % 3) plot the field U
 switch grid.d
  case 1
   for ie = 1:grid.ne
     my_pdesurf(h,Xi(:,:,ie),Ti,Ui(:,ie));
   end
  case 2
   np = size(Xi,2);
   nt = size(Ti,2);
   XXi = zeros(2,np*grid.ne);
   TTi = zeros(4,nt*grid.ne);
   UUi = zeros(np*grid.ne,1);
   for ie = 1:grid.ne
     XXi( : ,(ie-1)*np+1:ie*np) = Xi(:,:,ie);
     TTi(1:3,(ie-1)*nt+1:ie*nt) = Ti + (ie-1)*np;
     UUi(    (ie-1)*np+1:ie*np) = Ui(:,ie);
   end
   TTi(4,:) = 1;
   my_pdesurf(h,XXi,TTi,UUi);
  otherwise
   error('This functions supports only 1D and 2D.')
 endswitch
 axis equal

 hold off
 
return

